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Dual Scattering Channel Schemes Extending 
the Johns Algorithm 



Abstract. Dual scattering channel schemes extend the transmission line ma- 
trix numerical method (Johns' TLM algorithm) in two directions. For one point, 
transmission line links are replaced by abstract scattering channels in terms of 
paired distributions (characteristic impedances are thus neither needed, nor in 
general defined, e.g.). In the second place, non-trivial cell interface scattering is 
admitted during the connection cycle. Both extensions open a wide field of appli- 
cations beyond the range of classical time domain schemes, such as Yee's FDTD 
method and TLM. A DSC heat propagation [diffusion] scheme in non- orthogonal 
mesh, wherein heat sources are directly coupled to a lossy Maxwell field, illus- 
trates the approach. 
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Only a man who knows nothing of reason 

talks of reasoning without strong first principles. 

Gilbert Keith Chesterton 

1. Introduction 

It seems a paradox - and is just a typical process in mathematical analysis 
that a structure turns simple in a more general setting which at the same 
time widens its range of application. Accordingly - not very surprising - 
some ill-famed 'intricacies of the propagator approach to TLM' (sic. Rebel 
0, p. 5) virtually vanish if some of its elements are taken as the building 
blocks of a more general scheme. In fact, constructing the latter on essen- 
tially these elements in a quasi axiomatic manner will prove such intricacies 
to be mere artefacts of an inadequate framework. 

The choice of elements proposed in this paper 'generalizes' the Johns algo- 
rithm in two directions. In the first place, abstract scattering channels re- 
place transmission lines, which have some unpleasent properties (sectionEJ. 
Secondly, non-trivial cell boundary (interface) scattering is permitted dur- 
ing the connection cycle. The schemes thus obtained are characterized by a 
non trivial two-step (connection-reflection) cycle of iteration which exhibits 
certain duality relations - whence their name. 

When P.B. Johns and co-workers introduced the transmission line matrix 
(TLM) numerical method in the 1970s |2] it was almost instantaneously 
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assimilated by the microwave engineering community. In the same audience 
the method remained until today subject of assiduous study and extensive 
publication. Three conferences cxplicitely focussed on TLM [HElEl, and the 
monographs of Christopoulos [Sj and de Cogan [7] deal in detail with the 
original ideas as well as with classical applications. 

Familiarity with the transmission line picture, and the well known scat- 
tering concept, certainly fostered the acceptance of the TLM method among 
microwave engineers. On the other hand, just so the primary interest turned 
of course on applications in their own discipline, rather than onto the inner 
algorithmic structure as an object of mathematical analysis. Over the years 
still a node more was routinely invented, with new dispersion characteristics 
and/or equipped with still another ingenious stub, designed to model spe- 
cial propagation or transport phenomena, in varied geometries or boundary 
conditions. [HI stands somewhat exemplary for this line of research. 

Mathematical questions addressing the inner structure of the TLM al- 
gorithm and its potential generalizations have thus apparently been for a 
long time of secondary interest. They have yet not been left completely 
out of view. Chen, Ney and Hoefer [H] proved equivalence of the original 
(expanded) TLM node without stubs to the Yee finite-difference grid 

(TSUni. Recently, the non-trivial question of consistence of Johns' symmet- 
rical condensed node (SON), cf. Johns ^Jj, with Maxwell's equations and 
the intimately related problem of convergence to a smooth solution for de- 
creasing time step and grid spacing have been tackled, and in parts solved, 
by Rebel His thesis presents, by the way, a thorough survey over the 
ramifications of TLM until that time (year 2000), without perhaps spending 
sufficient attention to its non-orthogonal mesh extensions. 

From a quite general viewpoint, viz. widely independent of any par- 
ticular physical interpretation, the structure of the stub loaded (deflected) 
non-orthogonal TLM algorithm has been analysed in . The present paper 
goes even further and challenges the transmission line picture at all. The 
latter, in universally imposing free wave propagation between cells (with 
great benefit, at times), induces modeling limitations under circumstances 
that are outlined in section|2| Many restrictions can be by-passed by replac- 
ing transmission lines with abstract scattering channels in terms of 'paired' 
distributions. 

Dual scattering channel schemes are characterized by a two-step updat- 
ing cycle with certain duality relations between the two steps. The TLM 
method with its familiar connection-reflection cycle is trivial as a dual 
scheme in that the connection map reduces essentially to identity (viz. pure 
transmission or total reflection) - again with modeling limitations. These 
can be raised, anew, in permitting non trivial cell interface scattering dur- 
ing the connection step of iteration. 

One major merit of the transmission line method is unconditional sta- 
bility [12j . Since this property is usually drawn upon the passivity of linear 
transmission line networks, the question of stability needs a proper inves- 
tigation for DSC schemes that do not use hnes. That problem is studied 
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in |22j . where it is shown that a wide class of DSC schemes arc in fact 
unconditional stable. Due to the convolution type updating scheme {Johns' 
cycle; cf. section|21) it is sufficient for stability that the reflection and connec- 
tion maps share simple contraction properties (paraphrased as a-passivity 
in 1221). In summary, DSC schemes are unconditinally stable under quite 
general circumstances, and they are conceptually simple, though a set of 
technical definitions is of course ineluctable in a neat theory. Last but not 
least, nothing obscure nor 'intricate' should be associated anymore with the 
propagator approach. 

2. Scattering channels 

Any extension of the TLM method that includes heat transfer, fluid flow, or 
particle current, for instance, involves scattering channels other than trans- 
mission lines. The latter, for a non vanishing real part of the characteristic 
line impedance, inherently impose wave propagation between cells. Degener- 
ate lines, with a purely imaginary impedance, still work in diffusion models, 
cf. [H], chap. 7. Other types of transport or modes of propagation, such as 
for example the relativistic charged particle current treated in jl6j , are very 
unnaturally and more or less imperfectly modeled using transmission lines. 
There is good reason to get rid of lines in such and other cases within an 
extended framework. 

A first step towards the definition of more general scattering channels 
in TLM has been undertaken in replacing transmission lines with abstract 
projections into in- and outgoing field components, cf. |16| . It was postu- 
lated in this paper that the propagating fields ('link quantities') allow of a 
decomposition into a direct sum 

(1) Z = Zin ® Zouf , 

Zi„ and Zout representing the incident and outgoing fields, respectively. 
Moreover, it is essential in our understanding of TLM that the latter have 
a merely operational meaning in that only the total field z enters the dy- 
namical model equations (cf. sections I3I4|I . In singular physical 
interpretation can yet still be given to Zi„, Zout on the basis of a special 
analysis, e.g. ^Sj, Corollary 2. 

For the Maxwell field model the technical passage from the transmission 
line formulation to the projection operator setting is outlined in ^HJj Ap- 
pendix A. 

In the classical TLM setup the connection map simply transfers with- 
out further modifications the quantities outgoing from a cell into quantities 
incident at neigbouring cells or rejected at some totaly refiecting electric or 
magnetic wall. This is in perfect harmony with the behaviour of a propa- 
gating electromagnetic field, the components of which are tangential to the 
cell boundary (as the link quantities always are in a classical TLM cell, cf. 
PH) and that is thus not subject to refractive scattering at the cell face, 
even if the medium changes there. 
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The situation is clearly not thus simple for arbitrary propagating quan- 
tities. To circumvent any modeling restrictions, the connection cycle of a 
non-trivial DSC scheme comprises cell interface scattering from the outset. 
Nodal and cell face scattering thus enter a kind of duality relation that be- 
comes visible, for instance, in an apparent symmetry of the model equations 
in their most general form l|21l 1251 126() . Nodal and cell boundary scatter- 
ing may in fact be of equal importance and sometimes boundary scattering 
plays even the leading role in a DSC algorithm. 

In the generalized setup, just as in the traditional TLM framework, 
scattering channels interconnect a node, viz. a suitably defined centre of a 
mesh cell, with ports at the cell boundary. The channels are yet no longer 
represented by transmission lines. With respect to a computed physical field 
in D-dimensional configuration space, they simply form a pair of scalar or 
vector valued distributions, transposed over a distance in space, which test 
the field within the cell and on its boundary. A DSC scattering channel will 
thus be defined, precisely, as a pair of continuous linear functions {p, p"^ ) 
which act on a class of (suitably smooth real or complex) vector fields Z 
in configuration space, such that p has its support on a cell face and p"^ is 
connected to p via pull back into the node, i.e.: Given any notions of centre 
of cell and face, as well as the spatial translation s : MP that shifts 

the centre of a cell [node) into the centre of the face where p has its support, 
then the nodal image p~ of p is defined as the distribution 

(2) Z) := (po5, = (p, Zo,s-i) , 

and the pair (p, p~ ) is called a scattering channel. Equivalently, a scattering 
channel can of course be identified with ( p , s) or even simply with the 
port p, the pertinent shift and nodal image then being tacitly understood. 
The concept should in fact not be handled in too rigid a fashion - and 
there is no need to do so. In certain applications the support of the port 
distribution may be extended over a neighbourhood of a face, or the node 
distribution may rather be thought of as a mean over the entire cell (in 
the way familiar from finite volume methods). Needless to say that the 




Fig. 1. Ports on a cell face with their nodal images. 
stressed duality between nodal and cell boundary scattering is not to be 



Dual Scattering Channel Schemes 



5 



misunderstood in the narrow sense of category theory. Here, it refers simply 
to the observation that a set of propositions are vahd, modulo symmetry in 
certain terms, in the two scattering situations - which of course reflects the 
paired distribution concept of scattering channel and the already mentioned 
symmetry of the pertinent model equations in their most general form. A 
parallel symmetry then clearly characterizes the structure of the reflection 
and connection maps that solve these equations. 

Cell boundary scattering is. for the rest, not thus new an option: Already 
in the TLM model for superconducting boundary 201 cell face s-paramcters 
and boundary stubs have been introduced for solving the discretized London 
equations, cf. also [H^ . 

Despite the abolition of transmission lines, by their replacement with 
abstract scattering channels the computed ('physical') fields can still be 
represented - in the way familiar from the classical TLM method - as sums 
of in- and outgoing scalar or vector fields 

(3) Z = Ziji -\- ZquI 

No physical interpretation or propagation property is, however, in general 
ascribed to Zin^ out- In fact, these quantities are merely operationally defined 
by means of the well known Johns cycle of iteration 

Zin ( Q[Zout] + e ) 

(4) i t + T/2 

Ji and C denote the node and cell face propagators (or so-called reflection 
and connection maps - the latter including now cell boundary scattering, 
and e = e(t) induces any excitation. Note again that Zi„, Zout arc so far 
purely operational quantities, i.e. only the total fields z enter the model 
equations, while Zi„, Zout are in general bare of any physical meaning (a 
physical interpretation in terms of an energy flow still exists within the 
classical Maxwell field TLM model, cf. [Hj, CoroUary 2). 

As will be seen in the next section, the structures of the propagators 31 
and 6 are very similar in the general DSC scheme, thus reflecting the dual 
role that nodal and boundary scattering play therein. In a sense, precised 
in section 121 3? and C are the discrete convolution integrals that in every 
Johns cycle strictly solve the model equations. In fact, the Johns cycle can 
be looked at as basically a two-step convolution method for solving certain 
types of explicit finite difference equations in time. 

Note that the model equations can in principle be directly solved inas- 
much as they provide complete recurrence relations, cf. sections |21 0] The 
scattering approach (using Johns' cycle of convolutions) offers, however, 
important advantages. Thus, it provides clear cut reliable stability criteria 



T 

t + T 

T 



6 



Steffen Hein 



that make the DSC algorithm unconditionally stable under very general 
circumstances |22| . 

3. The elements characterizing DSC schemes 

So far, we dealt on a largely informal level with some typical traits of the 
TLM algorithm that either characterize DSC schemes in like manner, or 
which have to be modified in a specified way in order to attain a greater 
generality. We are still bound to keep our introductory promise and give a 
coherent description of DSC schemes in terms of some quasi axioms that 
condense their distinctive properties. Of course, we shall not really pursue 
axiomatics, here, in the sense of building a new theory on a complete set of 
basic assumptions. (Nor are we adopting a dogmatic attitude and going to 
fix a rigid framework that with certainty must be modified sooner or later in 
order to face a particular problem.) The emphasis is rather on compiling on 
a largely preliminary basis those formal elements that, in essence, lead to the 
peculiar structure of DSC schemes (in general), and of the TLM method (in 
particular), without being distracted by unnecessary information, such as 
mesh topology and geometry, s-parameters, etc., which characterize only a 
singular physical interpretation. A basic set of technical notions need simply 
to be clarified. 

In the following, 'simpZe' definitions are visualized in writing the defined 
object(s) in italics, more crucial or technical ones are explicitely designated 
as Definition. 

Until further notice, a mesh denotes a non-void finite family (i.e. an 
indexed set) of elements named cells ^ which are sets in their turn, and share 
the following properties. Each cell ^ contains an element , called its node, 
and a finite family dC = {dCi,}, called the (cell) boundary. The latter is built 
up of elements , named faces, which are sometimes simply written t in 
the place oi 8(1- 

Definition 1 A mesh M is called regular, if and only if it satisfies the 
following requirements of simplicity (S) and connectedness (C): 

(S) Every node belongs to exactly one cell and every face to at most two cells 
in M. 

(C) For every two cells CC'' G M, there exists a connecting sequence 

s = (C'')k=o G M^, such that = Ci C*^ — C"' o,^d, every two subsequent 
cells ^^^+1 in s have a common face, for < k < k. 

Also - certainly not too misleading: any two cells with a common face 
are called adjacent or neighbouring cells, and the common face a connecting 
face or interface. 

By a first postulate, DSC meshes are always regular meshes. 

The state space is any product of real or complex normed linear spaces 
labelled by mesh cells 



(5) 



s = TT §c 
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In addition, wc require that a DSC state space always contains a non-void 
subspace T C S , named the space of propagating fields , which on every 8^; 
reduces to a product of 'squared' spaces in the foUowing precise sense 



(6) ^c = 5'nsc = n,,,, 



In less formal language: every propagating field z splits over the cells into 
a sequence of pairs 

(7) = {z„ z';;'),(zoc 

labelled by the faces of the cell boundary. In other words, there is a canonical 
automorphism of normed linear spaces on J*, which on every J*^ reduces to 

nb : — > 

nb is obviously involutary ( — Id), and is called the node-boundary 
map. 

The components z^ and z"^ in Q. (jH)) are named the port (or face) 
component, and the node component, respectively, of z = (z^, 2~). They 
are usually written z^ = z^ = 7r^'(z) and = z^ = 7r"(z) with projections 
TT^ , tt" that are canonically extended over the entire space T. 

Let J : = UcGji/ be the set of all faces in M (remember that is 
defined as a union of faces). Then, in virtue of J' splits completely into 
subspaces 

(9) J'-n.,,^^ with V.:=U,.^^,,,P^X ■ 
A DSC process is a step function of time 

(10) pr:[0,T)^§ , 

such that TT^ o pr(t) and tt" o pr{t — t/2) are constant on every time interval 
[fiT, (fi + 1)t), /i g N , where they are defined. In other words, port compo- 
nents of a DSC process switch at integer multiples of the time step r while 
node quantities switch at odd integer multiples of t/2. 

Given a process, a state z with its entire history up to time t is usually 
written as a 'back in time running' sequence 

(11) [z]{t) (z(t-Mr))^gM , 

expanding so the domain of definition eventually to the negative time axis 
in the trivial way, i.e. z{s) : = for s < 0. 

By this convention, we assign thus to index /j, the (varying) state back 
in the past from present time t 

(12) [zUt) := z{t-fiT) , 
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rather than the (fixed) state z{fiT) - which has the technical advantage that 
fi so is directly related to a time difference (and eventually to an order of a 
finite difference equation in time), rather than to an absolute time (which 
is quite uninteresting, in general). 

Functions defined on back in time running sequences, such as are 
called causa/ functions (or propagators). Any such map is a discrete analogue 
to a causal Green's function integral, as has already been outlined in |18| . 

Let denote the set of all sequences with an arbitrary, but finite, 
number of non- vanishing elements in a linear space X. For every mesh 
cell C and face t in M consider then the subspaces of propagating fields 
'S>2 7r"(a'^), 7rP(J',). 

Definition 2 A refiection map (in £ M ) is a (possibly time dependent) 
causal operator 

and a connection map (in l £ J ) is a (likewise possibly time dependent) 
causal operator 

^ ' [^^ ^ e,[zP] . 

Also, a DSC system over M is a pair (6, 3?) consisting of any two families 

(15) e = {e,,},ej and 3? = {SlJ^eM 

of connection and reflection maps. Note that sometimes also the families 
are called the connection and reflection map. For the sake of algorithm 
stability, these maps (in whatever meaning) should share certain contraction 
properties studied in P^ . 

An excitation is only a distinguished process with values in the mesh 
boundary states. More precisely, let B : = {l \ l € J and t is not an interface} 
be defined as the mesh boundary, then an excitation is a process 

e : [0, T) — > TT 

(16) ^ ^^'^B 

t I — > e{t) ^ttP o e{t) , 

i.e. e is a port process, and hence switches at entire multiples of the time 
step T, and e generates non-interface {mesh boundary) states. 

Definition 3 The DSC process generated by (C, 31) and excited by e is 
the unique process z{t) = {z^, z"){t) which at every time t € [0,T) satisfies 

, . ^"(t) = nboz-{t + T/2) + zUt) 

^ ' z-{t) = z.m+nbozUt + rm , 
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the right-hand side being recursively defined through cyclic iteration of 



C(t + T/2) 

(18) Z^^,{t + T) 

t 



= n6o[e[2l,](t) + e(0 
- n6o3lK](< + r/2) 

= t + T 



(in that order) with, initially, ^/^(O) — z^^^{0) = 0. 
Remember that nb denotes the node-boundary map (|SJ). 

In lUHI) e and 3? stand of course for application of all propagators 
and 3^^ in the pertinent families (over J and M, respectively). Note that 
the order of application within the families is unimportant in virtue of the 
pairwise disjointness of all CPt, and of all T^- - which obviously implies that 
either C and are completely parallelizable at every time step. 

It follows immediately that z" and zj^ thus defined are node processes, 
hence switch at odd integer multiples of r/2, while z^, z^^^ are port processes 
(so they carry their superscripts aright). 
Equations H17|l . (|18|l can still be simplified to 

zv,n ^ ,P." + ,P;;^ and 

(19) zL{t) = Q[zLtm + m 

z:ut(t + T/2) = 3lK](f + r/2) , 

writing as usual zf^ and zJJ,j for 

zl{t) := n6oz,"„(t + T/2) 
n6ozP„,(0 =: zj;,(i-r/2) . 



(20) 



Comparing this to the TLM usage we note that in port quantities zf„, 
-^out ^'"c first introduced. With these are then node quantities zjjj, z"„( iden- 
tified (modulo the time shifts ±r/2, just as in H18|l ) without yet explicitely 
mentioning the node-boundary isomorphism nb. We shall sometimes follow 
this usage and omit the symbol nb where this cannot lead to confusion. 

Nothing has been said, so far, about physical interpretations or any 
implemented dynamical equations. In fact, the characteristic structure of 
the DSC algorithm is entirely laid down with the given definitions. As will 
be seen in the next section, typical features and facts, some quite familiar 
from TLM, are derived straight away with only the above elements. 

With respect to the dynamical model equations - which the algorithm has 
ultimately to solve and that determine the propagators 3?,^, 6^, cf. section0]- 
we reiterate the important general agreement that only total fields z^, z", 
not, however, their incident and outgoing components separately , shall enter 
these equations. Accordingly, we consider only model equations between 
total fields. Quite generally, and modulo further restrictions (inferred in the 
next section), the DSC model equations should be of the types 

(21) + 
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with causal functions and shortly [z±] := [z ]( i ± r/2 ). The r/2 

time shifts synchronize node and cell boundary switching in H21(l , such that 
the equations can be strictly solved, for every time t S [ , T ) and are well- 
posed, in this sense. Of course, time shifts by — t/2 would also lead to 
synchronization. The resulting equations would, however, conflict with the 
causality property of Jl and C. 

Inspection of equations (|21(l shows that 5"" affects only the reflection 
cycle, while 3^^ has impact only on the connection cycle. We are now dealing 
with the model equations in some more detail. 

4. The model equations 

The physical interpretation of a DSC system fixes, intuitively speaking, 
the terms in that states in T are read as physical fields. More deliberately, 
certain states in CP are interpreted as distributional values (finite integrals, 
e.g.) of physical fields, which are localized in a mesh cell system. 

Any interpretation requires, hence, in the first instance a geometric re- 
alization of the underlying regular mesh, wherein the relations between ab- 
stract cells, nodes, and boundary faces which characterize M are translated 
into relations between geometric objects, bounded subsets of R^, such as 
(in general) polyhedral mesh cells with their faces, e.g. 

Given any geometric realization of M, a physical interpretation of a DSC 
system over M is, precisely, a family / = ^gg^ of continuous linear 

functions 

(22) , 

each defined on a space of smooth m-component vector fields in D- 
dimensional configuration space (i.e. C C^{M.^)™; m € N depending on 
(i, C) ), such that has its distributional support on a cell face and its 
range in CP^ = tt? (?<;). 

Note that index l in H22(l may optionally be read as a cell face label or 
as a multiindex referring to a set of ports on the same face. Since we are 
dealing with vector- valued distributions (with range in T ) and the support 
of every is required to be localized on a cell face (which can be weakened 
to at least associated to a face) , there is essentially no difference in reading 
zj'^ — I'q (Z) as a cell face state vector, or as an array of components (labelled 
by port indices) of such a vector. Thus, index l in (|22|l may be thought of 
as implicitely labelling a subset of ports on face l E d<^. 

Attention is also drawn to the fact that the functions are not required 
to be surjective onto 7^ (i.e. not every state in must be directly related 
do a distribution in /). There is, for instance, no need to exclude from 
any function or linear combination of fields in different spaces £' (which 
may represent a spatial finite difference of fields, as in the approximate 
gradient of our sample model in section |S1). 

The evaluation of nodal fields goes quasi pick-a-pack with / by applying 
the scattering channel concept of sectional 
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Let sj^ denote the translational shift in R-'^ from any node into the 
(centre of the) face where has its support, of. figG] and assume without 
loss of gencrahty that : Z{x) i— > Z{x+s''q) is an inner map in (otherwise 
take the closure of under such transformations). Then clearly holds 

Proposition 1 For every I'^ E I there exists exactly one function I'^ : — s- 
( note Iq ^ I ), such that the following diagram is commutative 



nh 



(23) I I If 



£' > £^ 

i.e. If oS'^=nboI^. 

Proof Mere retrospection of definitions. 

In the terminology of section [21 is the nodal image of the port(s) I^, 
and the pair of distributions { If nb o ) forms a scattering channel. 

The dynamical DSC model equations are, in the line of the preceeding, 
to be read as finite difference equations in time between states 
z = izf,zf) G y that have an interpretation as distributional values of 
physical fields Z G £J. , 

(24) <c = ^c(^) ' <C=n(Z) . 

Unlike classical FD equations between pointwise evaluated physical fields, 
the DSC model equations interrelate in many cases finite integrals over 
a line segment or face, e.g. - pointwise evaluation (with a Dirac measure 
as distribution) not excluded. The DSC approach is, in this respect, by 
far more versatile than the classical finite difference time domain method. 
So, the former permits, for instance, of generalizing to a non-orthogonal 
mesh Johns' TLM method JT] in much a simpler way, cf. ^Tj, than the 
FDTD approach allows for Yee's method of approximation to Maxwell's 
equations 

A central principle underlying DSC schemes is near-field interaction. 
Like causality, this is already implicit in the (domains of) definition of the 
reflection and connection maps. Near-field interaction simply spells that 
only the fields in the immediate neighbourhood of a state z - precisely only 
those in if z G Tf and those in T,, if z G , cf. P EJ, along with 
their history, of course - determine the evolution of that state on the next 
updating step (note, this refers to fields evaluated in 7' and not, for instance, 
to an exterior potential, which may still induce a time dependence of Ji or 

e). 

In other words, an updated nodal state depends only on states of the 
pertinent cell, including its boundary, while the evolution of a port state 
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is determined by states of the respective face and by nodal states of the 
adjacent cehs. 

It follows that the model equations (|21|) split into the two families 

(25) ^cK+H^] ^ 0, z^e^c, CeM 

(26) ^nz^+]K] ^ 0, 2,, eT,, , LeJ. 

(Remaining aware of the dependence upon C and t of these equations, we 
can in general omit the subscripts, if there is no danger of confusion.) 

Since the following analysis runs perfectly parallel for the dual equations, 
from now on it is confined to the implications of H25|l . - The reader may 
write down parallel statements for dual equations, at times, by exchanging 
port for nodal and incident for outgoing quantities, starting, for instance, 
with the cell boundary version (with C in the place of 31) of the following: 

Definition 4 Let I ^ [Q, T) be a finite intervall (i.e. T G M+ ). Then we 
shall say that "R generates solutions of the model equations (|25() on /, if 
and only if for every sequence of incident nodal fields [z"^] the (obviously 
unique) process z = Zi„ + Zout which is recursively given by 

zlit - r/2) = 

(27) z:,,{t) = 
4ut{t + r/2) = 

algebraically solves equations H25(l (identically on I). Sometimes, we then 
simply say that 5i solves the model equations on that ( finite !) interval. 

Remark 1 

(i) Note that Definition 01 refers to a purely algebraic property of 3^ that 
is not yet related to any stability questions, for instance (cf. also the 
Remark to Theorem n). 

(ii) If 3? solves equations (|25|l on /, then in particular every process generated 
by (6, 3?) in the sense of Definition 0] solves H25|l on /, since every such 
process obviously satisfies H27() (ef. (jSOl)- 

The least awkward (yet fortunately frequently encountered) situation is 
brought abount with homogeneous linear equations, i.e. for the evolution of 
nodal states 

(28) ?"[z^][2f] =^J^^0^z"(t + T/2-^r)+^^zf(t-^T) = 0, 
with linear and possibly time dependent operators 



nbzlXt) 
nbz-{t) 
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that map into any common linear space 3 wherein 3^ has its range. Similar 
dual equations 5'P[z?jl] [z"] = determine the linear evolution of the port 
states during the connection cycle. 

Many, if not almost all (viz. all but a finite number) of the 0^, may be 
zero. Any maximum /.j e N U {oo}, such that 0^ 7^ or -0^ 7^ 0, is called the 
(dynamical) order of the model equations, and in general equals the order 
in time of the integro-differential equations which physically describe the 
underlying dynamical problem in terms of smooth fields. 

The following statement provides a theoretical means for computing the 
reflection map of equation (|28|l recursively. 

Theorem 1 Let (j)Q : — > 3 be bijective (i.e. one-to-one and onto). Then 
solves equation (|28(l on a finite interval [0, T), if and only if for every 
tG[0, T) 

mjit) = z:^,{t) = 

+ (0/.+1 + i^t^nb) z^^^it -T - fir) } . 

Proof Substituting in (|28|l for z P' " the right-hand sides of H17|l and using 
equations H2(J|I yields for t < T on a time domain trivially extended over 
the negative axis 

{ '^'^ [ ^™ (t + r/2 - ^1T) + (t + t/2 - mt) ] + 

+ ij^nb[ zZ, {t + T/2- ^lT) + zl, (t - r/2 - /ir) ] } ^ , 

which for invcrtible 0o and t' : = < + t/2 is equivalent to the recursion for- 
mula of the theorem. 

Remark 2 For every t < T < 00 the sum in Theorem ^ is actually finite, 
hence convergence is not a question as long as one abstains from considering 
limits. For finite T , the theorem conveys purely algebraic relations inherited 
from the structure of the DSC process as laid down by (|17II18|1 in section|3 
The question of stability of a process for T — » cx) is being addressed below. 

If 00 is not bijective, then the model equations H28|l are incomplete in 
that they do not determine a well-defined and unique reflection propagator 
in every cell. We therefore require that the uniqueness conditition 

(U) 00 : 5*" ^ 3 is bijective ( i.e. a linear isomorphism ) 
be always satisfied. Clearly, l)28|l then defines an explicit scheme. 

The model equations of the classical TLM method discretize Maxwell's 
equations. So, they are linear and first order in time, viz. of the type 

(29) 00 z"(t + r/2) + 01 z"(t - t/2) + 0o z^it) , 

with time independent operators 0o. 1 and -00, which are derived in |17l ?]. 
for instance. Of the same type - viz. linear and of first dynamical order - 
are the discretized diffusion equations. Theorem ^ applied to this special 
situation yields 



14 



Steffen Hein 



Corollary 41 The (for a given time step unique) reflection map that solves 
on a finite interval the first order linear equations H29() with time indepen- 
dent (/) 0, 1 and ipo is 



(30) 

wherein 



K = —Id — (j)^^ ijjQ nb 
L = -00-1 

M = + 01 + ( 01 + Vo ) ( - 00 M( <?^o + V'o '^^ ) 



K 

Proof By induction, applying Theorem^ to an incident Dirac pulse 

^Ut) ■■ = ZX[0.r)(t-T/2) 

with any fixed state vector z G IP" ; Xi denotes the characteristic function 
of interval / (which equals 1 for every argument in / and elsewhere). By 
linearity the statement then holds for arbitrary incident processes zf^{t), 
each of them being writable as a superposition of Dirac processes that start 
at subsequent time steps. 

Remark 3 A sufficient condition for convergence of the propagator series (|3U|I 
(applied to any finite incident pulse zf^ ) in the limit t — > cxd is obviously 

(31) II II < 1 , 

where || ... || is any submultiplicative norm, e.g. the Hilbert {spectral) 
norm 

(32) II N \\h : = max { | A | ; A eigenvalue of N }. 

In fact, the latter condition is sufficient for algorithm stability of the Maxwell 
field TLM model, as shown in |21II17| . Any first order linear process is clearly 
stable, if the Hilbert norms of K, L, M, N are bounded by 1 (strictly for 
K and N), since the propagator 31 then is contractive. This is in general 
ensured with bounds for the time step (if necessary, in combination with a 
transformation (|34|l '). 

For linear model equations of any finite order, recursion formulae that 
generalize (|30|) are easily derived from Theorem ^ e.g. ^Hli equation (32). 

Corollary 42 With K, L, M, N as above, the DSC process solving 
permits a representation as a 'deflected' scattering process 

Ct(0\ ^ fK L\ ( z^^{t) 
d{t) ) \MN) \d{t-^iT) 

the deflection d{t) being completely recursively defined with the (universally 
maintained) initial conditions rf|t<o = 0- 



(33) 
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Proof Immediate consequence of Corollary 1411 

Remark 4 Note that there remains an arbitrariness in the definition of the 
operators L, M , and N in Corollaries l41l In fact, any invertible trans- 
formation 1:3—^3 together with the simultaneous replacements of L, M , 
and N by, respectively, 

(34) L~:=Lo/-i, M'" := /o7\/, and : = / o TV o /-^ 
clearly do not alter the propagator 31, i.e. generates the same DSC process. 

The corollaries offer solutions of first order linear equations with time 
independent operator coefficients which are complete in the sense of (U). 
They thus cover the entire field of classical TLM (with connection maps 
that are trivial in reducing essentially to identity). 

In certain situations it may be useful, or necessary, to integrate some 
new (possibly non-linear) interactions into a given DSC model. Sometimes 
this can be carried out by adding suitable coupling terms to the equations 
of the yet existing model. We therefore consider perturbed model equations 
of the type 

(35) y~[zlW] = 7''[zl][zP]+3[z''_][zP] = 

(here exemplary for nodal perturbations), wherein 3 denotes any causal map 
into 3. The (— r/2) time shift in the first argument of 3 again synchronizes 
port and node switching. Note that the time shift is negative, here. This is 
to ensure that the perturbation 3 cannot destroy the uniqueness conditions 
(U) in the case of a linear function 3^", and to preserve expliciteness of the 
updating relations, in general. The importance of this condition becomes 
clear in the proof of the following formula. 

Proposition 2 (Deflection Formula) Let the reflection map generate 
solutions of the linear equations H28|l on a finite interval I — [0, T). Then 
IK solves the perturbed equations H35() on I , if and only if the so-called 
deflection 

X) := 3l~ - 3? 

satisfies recursively 

Proof By definition, 3i solves equations (|35ll on /, if and only if for every 
incident sequence [ ] and 

Zout ■ ~ \Zin\ J Z : = Zin + Zout 

holds 

^^[zi][zn = 

= J"[zl][zP] + 3[z1][zP] = . 
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However, with 

Wont ■ ^[Z?n\ : W : = Zi^ + Wont 

( hence T> = z^^^ — li;"^^ ) and using (^5)1 , this is the case iff 

1^^'-^ V ^ ■' 

+ V'a.I* (^I'„|t+r/2-,.r + K,ut\t-rl2-i,r + ^3 1 t-r ) } ■ 

^OUt|t-T/2-,.lT 

In virtue of the linearity of the 0^ , , the latter identity holds iff 



by 



+ 0O|t 2^|t+T/2 + JZ^eN ('^'^ + 11* ^^1* ) 2^|t-r/2-/iT , 

which is the recurrence relations of the proposition. 



Corollary 1 (Deflected processes) Let 31 solve (|28|l on a finite interval 
I and (po '■ y "J be any bijective operator {that thus satisfies the complete- 
ness conditions (U)) . Then 

25|t+r/2 -.^ ^^-^]{3[Z1][Z% + Y.^^n^^^+'\' + V'M|t)25|t-r/2-^r } 

with initial conditions 'D|t<o = defines recursively a causal operator D, 
such that "R : = + D solves equations H35(l on I . 

Concluding this section, we stress once again that Theorem ^ and the 
ensuing propositions and corollaries apply just as well to the connection 
cycle, i.e. to cell interface scattering, provided the replacements (j25|) by 
3? by C, by Tt, port by node superscripts, and incoming by outgoing 
fields (and vice-versa) are simultaneously made. - Note, however, that any 
excitations may temporarily violate the model equations at a mesh boundary 
face. The model developer is encouraged to care for physically consistent 
excitations. 



5. A heat propagation scheme in non-orthogonal mesh 

The physical interpretation underlying the following application relates a 
smoothly varying (viz. in time and space continuously differentiable, C^-) 
temperature field T, evaluated as at the face centre points and as T" 
in the nodes of a of non-orthogonal hexahedral mesh, to total states 2^'" 
of a DSC model. In fact, we derive the model equations for the connection 
and refiection cycles of a DSC heat propagation (diffusion) scheme. Since 
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the equations are linear and of dynamical order and 1, respectively - as 
will be seen - they can be processed, following the guidelines of the last 
section. In the end, we display some computational results of a dispersion 
test carried out with this model. 

In order to simplify the notation we follow Einstein's convention to sum 
up over identical right-hand (!) sub and superscripts within all terms where 
such are present (summation is not carried out over any index that also 
appears somewhere at the left-hand side of a pertinent symbol - thus, in 
{—l)'^ a'^bx i^c the sum is made over A but not over k). The shape of a 



e 




re (a) (b) 

Fig. 2. Non-orthogonal hexabedral mesh cell. 
(a) Edge vectors. (b) Node vectors. 
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Fig. 3. Face vectors and temperature points (nodal section). 



hexahedral cell is completely determined by its 12 edge vectors (i/e)^=o....,ii- 
Also, with the labelling scheme of figE^, node vectors (;^&)^=o,i.2 and face 
vectors (t/)t=o,...,5 are defined as 

1 

ub := - > ,^ ^ ,e A* = 0,1,2 

^ ' and J := ( (8+.„ ^ ^ ) A 

A ((4+20^ +(5+20 e) t = 0,...,5 
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(edge vector indices cyclic modulo 12 and A denoting the wedge ('cross') 
product in M'^). 

At every face ie{0,...,5}ofa mesh cell and for any given r e M+ 
the following time shifted finite temperature differences in directions ^5 
( /i = 0, 1, 2 ) form a vector valued function 



(37) 

r 2 (-!)'( T" I ,_,/2 - ,Tf I , ) if = [l/2] 



.V^T^(t,r; 



{[x] denotes the integer part of a; G M ) . The time increments are chosen to 
attain technical consistence with the updating conventions of DSC schemes. 
They do not destroy convergence, as easily seen: In fact, in the centre point 
of face L the vector tV^T approximates in the first order of the time 
increment r , and of the linear cell extension, the scalar products of the 
node vectors with the temperature gradient VT . Let, precisely, for a fixed 
centre point on face t and e G the e-scaled cell have edge vectors 
tC'^ : = e tC . Let also ^V^ denote function (|37|1 for the e-scalcd cell 
(with node vectors ^6"" = ). Then at the fixed point holds 

(38) < ^6,grad(T) > = ■ VT = lim lim i.V^^T^, 

as immediately follows from the required C^-smoothness of the temperature 
field T. 

To recover, in the same sense and order of approximation, the gradient 

VT from H37() . observe that for every orthonormal basis {^u)^^^) „i-i of 

V = M™ or C™ , and for an arbitrary basis (^6)^=o,...,m-i with coordinate 
matrix = < , ^5 >, the scalar products of any vector a ^ V with ^6 
are 

(39) < ^b,a > = 

= (^2) = (/3?) 

(At the right-hand side, and henceforth, we follow Einstein's convention to 
sum up over identical sub and superscripts within terms where such are 
present). It follows that 

(40) a. =7>^ with {lt^):^m)T' ■ 

Loosely speaking, the scalar products of any vector with the basis vec- 
tors ^& transform into the coordinates of that vector with respect to an 
orthonormal basis by multiplication with matrix 7 ~ (/3*)^^ , where 
= < , p6 > , i.e. [3 is the matrix of the coordinate (column) vec- 
tors fj_b with respect to the given ON-basis , and 7 is the adjoint inverse 
of p. 
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This applied to the node vector basis and H^{8|l yields the approximate 
temperature gradient at face l as 

(41) = Jti .V^T^. 

Thus, the heat current into the cell through face i with face vector 
components i,/" = < ^.f , i^u > , e {0, 1, 2} , is 

Xh denoting the heat conductivity in the cell. 

The heat current through every interface is conserved, i.e. between any 
two adjacent cells C, x with the common face labelled l in cell and k in x 
applies 

(43) U = -IJ. 

Also, the nodal temperature change in cell C is 

where Cy denotes the heat capacity (per volume), V the cell volume, and S 
any heat source (s) in the cell. 

We finally introduce quantities tzP'" (t = 0, ...,5; /i = 0, 1, 2 ), which still 
smoothly vary in time with the temperature T (and that are hence not yet 
DSC states, but will later be updated as such) 



(45) a^it) 
and 

(46) ,zP{t) 



2(-i)'r"|t i^^i = [L/2] 

TP- 2f.TP)\t-r/2 else 



2(-l)Srf|, 

,,z;'(t - r/2) else. 



From (123 EH) follows 

.J I t+r/2 = ( I * - 2 (-1)'<5[:/21 ,TP I ) 



(47) 



This, with (|43|l and continuity of the temperature at the interface, ^T^ 
iTP , together with ^ imply 
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which form a complete set of recurrence relations for (given z") and 
so can be taken as model equations for the connection cycle of a DSC 
algorithm. 

Equations (|44|l discretely integrated in the time balanced form with 
increment r yield 

(49) = r"|t-./2+^(^+ ELo'-^I*)' 

i.e., with (|45l 1461 ITzjl . the recurrence relations 



else , 

which provide a complete set of model equations for the reflection cycle of a 
DSC algorithm. Note that the first line, modulo the factor 2 (—1)'' , always 
updates the nodal temperature. Of course, this has to be carried out only 
once per cell and iteration cycle. In this - typical - example the dual state 
space concept of DSC ( needed by Johns' cycle ) creates a redundancy, which 
is a price for process parallizability within either parts of the connection- 
reflection cycle. 

Equations gSl EOI) 

can be directly taken as updating relations for total 
quantities of a DSC scheme. Alternatively, they may be further processed in 
deriving reflection and connection maps, along with stability bounds for the 
time step. The proceeding is canonical and amounts in essence to a straight- 
forward transcription of the model equations along the lines of Theorem ^ 
and corollaries. 

It is particularly easy to couple this heat conduction model - within one 
and the same mesh - to a Maxwell field TLM model in the non-orthogonal 
setting In fact, with the node vector definition in ^HIj the total node 
voltages are just the scalar products of ^6 with the electric field, hence the 
dielectric losses and heat sources per cell are 

(51) S = \aVE'^E: = \<rVY.J^iU-W 

for a frequency domain (complex) TLM algorithm, cf. ; cr = 27r/ e tan((5) 
denotes the effective loss current conductivity at frequency / in a mesh cell 
of absolute permittivity e and dielectric loss factor tan((5) ; 7 = (/3*)~^ as 
in H40() . In Spinner's Maxwell field solver the model couples, in addition, 
to magnetic and skin effect losses. 

Finally, Fig 01 displays the result of a dispersion test, computed in a 
square mesh using non-orthogonal cells. It turns out that the heat conduc- 
tion properties of the mesh arc highly insensitive to cell shape and orienta- 
tion (as of course should be the case). 



(50) 



\ t+T/2 



[t/2] I t"T/2 + 2c„y 1 "5 + 
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Fig. 4. Transverse heat conduction over a square mesh using non-orthogonal cells. 
A Heaviside temperature step is imposed on one side and the transient tempera- 
ture computed at the opposite side, assuming adiabatic boundary conditions on 
all but the heated sides. 

DSC results (+) are plotted over analytical solution (curve). 
(a) The mesh, (b) Horizontal (c) vertical propagation. 



6. Conclusions 

Johns' TLM algorithm can be extended with benefit in two major directions 
by replacing transmission line links between cells with abstract scattering 
channels in terms of paired distributions and in admitting non-trivial cell 
interface scattering. Executing this program lead us in this paper to a new 
class of Dual Scattering Channel schemes which offer enhanced modeling 
potentiality and canonical techniques for stable algorithm design. 
spinner's implementation of a heat propagation scheme coupled to a 
lossy Maxwell field illustrates the approach. 

The connection and reflection cycles of a DSC process are (either) com- 
pletely parallelizable, which can be turned into account in computational 
performance. DSC schemes open a challenging field to future research. Ap- 
plications to fluid dynamics are presently under examination. 
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